﻿#light

open dnAnalytics.FSharp
open dnAnalytics.LinearAlgebra
open dnAnalytics.Statistics.Distributions

/// The mean vector of the distribution.
let m = vector [2.5; -5.0]
/// The covariance matrix of the distribution.
let C = matrix [[1.0; 0.9] ; [0.9;1.0]]

/// The random number generator to use.
let rnd = new System.Random()
/// The Cholesky factor of the covariance matrix.
let L = (new Decomposition.Cholesky(C)).Factor()
/// The multivariate samples.
let samples =
  Array.init 100  
    (fun _ -> vector [Normal.Sample(rnd, 0.0, 1.0);  
                      Normal.Sample(rnd, 0.0, 1.0)])  
  |> Array.map (fun v -> L * v + m)


/// The empirical mean of the samples.
let eMean =
  samples 
  |> Array.fold_left (+) (vector [0.0;0.0]) 
  |> (fun x -> x / 100.0)

/// The empirical covariance of the samples.
let eCoVar =
  samples 
  |> Array.fold_left 
      (fun acc v -> acc + (v - m) * (v - m) / 100.0) 
      (matrix [[0.0;0.0];[0.0;0.0]])


printfn "Empirical mean is [%f,%f]; true mean is [%f,%f]" eMean.[0] eMean.[1] m.[0] m.[1]
printfn "Empirical covariance is [[%f,%f];[%f,%f]]; true covariance is [[%f,%f];[%f,%f]]"
    eCoVar.[0,0] eCoVar.[0,1] eCoVar.[1,0] eCoVar.[1,1]
    C.[0,0] C.[0,1] C.[1,0] C.[1,1]